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Abstract 



Q ! A cluster Monte Carlo algorithm for the Ashkin-Teller (AT) model is con- 

I structed according to the guidelines of a general scheme for such algorithms. 

■<!::;j- ' Its dynamical behaviour is tested for the square lattice AT model. We per- 

form simulations on the line of critical points along which the exponents vary 
continuously, and find that critical slowing down is significantly reduced. We 
, find continuous variation of the dynamical exponent z along the line, follow- 

' ing the variation of the ratio a/u, in a manner which satisfies the Li-Sokal 

' bound Zciuster ^ Oi/v, that was so far proved only for Potts models. 
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I. INTRODUCTION 



The Ashkin- Teller (AT) mo del has been studied extensively, ever since its introduction, 
by a variety of methods. In two dimensions, in particular, much is known about the phase 
diagram and critical behavior of the model . Nevertheless there are problems that were 
not addressed extensively before; the critical behavior of random AT models is one such 
issue about which not much is known, hope for analytic treatment is slim, and therefore one 
expects numerical simulations to be the main tool of investigation. With this aim in mind, 
we set out to develop an efficient Monte Carlo (MC) cluster algorithm for the AT model. 

A convenient representation of the generalized AT model (gAT) is in terms of two Ising 
spin variables, cij and Tj, placed on every site of a lattice. Denoting by < ij > a. pair of 
nearest neighbor sites, the Hamiltonian is given by 

<ij> 

Here Kcr {Kr) are the strengths of the interactions between neighboring a (r) spins, and L 
is a four-spin coupling. The phase diagram of the ferromagnetic general AT model is known 
in two dimensions from duality transformations and renormalization group studiesili. The 
three dimensional model has been studied as well@. In this paper we are concerned with 
the Z{4) subspace of the general model, in which Ka^ = Kr = K. The phase diagram in 
this subspace is reviewed in Sec. 2. The critical properties of the model in this subspace 
are of special interest, since it has a line of critical points, along which the exponents vary 
continuously, and have been determined analyticallyO, interpolating between Ising and four- 
state Potts exponents. For instance, the value of the ratio a/v varies from at the Ising 
(L = 0) critical point to a/i/ = 1 at the four-state Potts point K = L. This exponent is 
of special interest to us since it has been proved that for Potts models it serves as lower 
bound to the dynamic exponent z of cluster algorithm^ Furthermore, the bound seems to 
be reached in the case of the Ising and Potts mo delsi. 

Cluster algorithms0 as introduced by Swendsen and Wangi, and extended by Wolfi, are 
reviewed briefly in Sec. 3. These algorithms give rise to dynamics whose relaxation time, 
Tswi is significantly shorter than that of standard, single-spin flip MC methods. This is 
most important at a critical point, where the relaxation time of a finite system grows with 
its linear size L according to 

r ~ . 

Cluster algorithms have a significantly lower dynamic exponent zsw than that of standard 
MC. Hence if one is interested in performing extensive simulations of a model such as AT, 
it is well worth to spend time on developing an appropriate cluster algorithm. 

Creating an efficient algorithm can be a challenge. Naive application of the original SW 
scheme does not work (as explained in Sec. 3 and demonstrated in Sec. 4). We set out 
to construct an efficient cluster algorithm using the guidelines and methods that were pre- 
sented by Kandel and Domany0 . This general scheme is guaranteed to yield an algorithm 
that satisfies detailed balance, and once the important excitations of the model have been 
identified and incorporated, we are guaranteed to get a good algorithm. This general for- 
malism is briefiy reviewed, and the resulting algorithm is presented in Sec. 3. Interestingly, 
the cluster algorithm found this way is identical to what would have been obtained had we 
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used Wolff's Ising embedding methodli0, as shown in the Appendix. Numerical results are 
given in Sec. 4; in particular, efficiency of our algorithm is tested by measuring the dynamic 
exponent along the critical line, which is indeed significantly lower than that of standard 
MC. 

An interesting question we set out to resolve concerns comparison of zsw with a/u along 
the AT critical line. We found that the Li-Sokal boundi 

a 

^cluster ^ 

V 

is satisfied, variation of the dynamic exponent follows that of a/i^, and within our numerical 
accuracy and limitations due to finite size effects, our results indicate that the two are equal. 



II. PHASE DIAGRAM OF THE ASHKIN TELLER MODEL 

We now review the phase diagram of the square lattice AT modeli and some of its 
critical properties. Figure |1| gives the phase diagram of the model plotted as a function of 
the parameters X and Z where 

Z = exp-4^ ; X = exp-2(^+^) . (1) 

At times we use the notation X = {X, Z) to denote a point in phase space. Full lines of 
phase transitions separate three phases: 

(a) Paramagnetic, labeled as "P" . The couplings are sufficiently weak so that the system 
is in a paramagnetic phase in which neither a nor r (nor ar) are ordered. 

(b) Ferromagnetic phase, labeled as "F" . The couplings are sufficiently strong so that 
cr and r independently order in a ferromagnetic fashion so that (cr) = ±(t). In this phase 
(ar) is also different from zero and has the same sign as (cT)(r). 

(c) A phase labeled (err), in which cxr is ordered ferromagneticaly but (a) = (r) = 0. 
This phase arises only for L > K. 

On the dashed line Z = X^ we have L = ; and thus it is obviously a subspace of two 
decoupled Ising models, having an Ising transition at the point Xq. The dashed line Z = X 
has L = K, in which case the AT model becomes the 4-state Potts model. The point X4 is 
a 4-state Potts multicritical point. 

A marginal operator generates a continuous variation of critical exponents along the 
line X0X4, isomorphic to the known critical line of the eight-vertex model0. Through an 
exact duality type transformation this line is mapped onto the critical line of a staggered 
8-vertex model and through a relation with the Coulomb gas its critical exponents are known 
exactly!!: 

yt = 2- {2/gR) (2) 

where 

(7R = -sin-i[icoth(2i^)], (3) 

and ^ = I all along the line. Lastly the lines X^B and X4C flow under renormalization to 
Ising type fixed points. 

The exact location of the transition line XqX^ can be found through the duality trans- 
formation of the AT modeli. It is given by the self-dual line Z = \ — 2X. 
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III. CLUSTER METHOD FOR THE ASHKIN TELLER MODEL 



A. SW cluster method 

Cluster algorithms have proved to be a useful method of reducing critical slowing down 
in MC simulations. For completeness we review here the pioneering cluster algorithm of 
Swendsen and Wang (SW)i for the Ising modelEl with the Hamiltonian 

n = -j ^fc^i (X = ±1 . (4) 

<k,j> 

The SW procedure stochastically identifies clusters of aligned spins, and then flips whole 
clusters simultaneously. Starting from a given configuration u, SW go over all the bonds, 
and either "freeze" or "delete" them. A bond connecting two neighboring sites k and j, is 
deleted with probability P^, and frozen with probability Pf = 1 — P^, where: 

= e~P{J-,-,+J) . (5) 

Having gone over all the bonds, all spins which have a path of frozen bonds connecting 
them, are identified as being in the same cluster. Now the new configuration is generated 
by fiipping every cluster with probability 1/2 . Note that according to (^), only spins of 
the same sign can be frozen in the same cluster. SW identify correctly the elementary large 
scale excitations as clusters of aligned spins. This correct identification is essential to the 
success of a cluster method for any more general model. 

In sec. |1IID| we shall present a naive implementation of the original SW scheme to the 



AT model, and explain why it is not expected to work well. This necessitates search for 
a different clustering rule; the one we found is based on a general scheme which we now 
describe. 



B. General scheme for cluster methods 

A unifying view of all Cluster Algorithms has been given by Kandel and Domany (KD)0, 
which we now review. The general scheme consists of two steps: given a spin configuration 
u, the first step consists of stochastically generating a new Hamiltonian 7i. The second step 
consists of simulating the model with the new Hamiltonian, thereby bringing it to a new 
configuration u'. To carry out the first step, KD write the Hamiltonian as 

n = T.Vi- (6) 

I 

For example Vi can be the energy of a single bond in the nearest neighbor Ising Hamiltonian. 
Then to each they assign one of possible integers i, in a stochastic manner that depends 
on the starting configuration u. That is, the probability to assign z to / is written as 
Pi = Pliu). This probability is normalized, i.e. 

E^'H = 1 (7) 
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for any term / and configuration u. Tlien they construct a new Hamiltonian: 



^W = EK(0, (8) 

I 

where (for any spin configuration u): 

Vl{u) = V,{u)--^HPl{u)]+C\. (9) 

The free parameters C\ are configuration independent. 

The second step consists of simulation of the model with any procedure whose transition 
probability T{j}(M — »• u') satisfies the detailed balance condition with respect to the new 
Hamiltonian, i.e. 

e-/3^o}Wf|,.j(y _ u') = e"'^^{'>("')f{i|(M' u). (10) 

After completing the second step, a new configuration u' is arrived at, the original Hamilto- 
nian is restored and the process is repeated. Equations ( and ( |10D ensure that the whole 
procedure satisfies the detailed balance condition, with respect to the original Hamiltonian 
(for the proof se ), but ergodicity needs to be proved for each application separately. 

We give now two types of modifications to the Hamiltonian. Consider a term Vi{u) that 
can take M distinct energy values Ei, i = 1,. . . ,M. The first is the deletion operation, used 
by SW, which eliminates the interaction Vi{u) that gets replaced by 

yj(n)=0, (11) 

for any configuration u. To get this, we must have ( see eq. ( 

pj(«) = e^W(")+^^]. (12) 

The constant must be chosen so that -?](«) < 1 for any u. 

The second modification is a 'generalized' freezing operation which we will later use in 
our scheme for the Ashkin- Teller model. Its probability is 

pi(. ^ f e^W(")+^^] if Vi{u) = E„ and 1 < m < /i . . 

^ \ if Vi{u) = E„ and /i < m < M ' ^ ^ 

and the modified interaction according to (^) is 

VU~\ = I^ if V^(m) = and 1 < m < /i , . 

9^'^' \ oo if Vi{u) = and m < M ■ 

Again the constant Cg must be chosen so that Pg{u) < 1 for any u. This type of modification 
allows free movements between some of the states u of Vi{u) which didn't have the same 
energy in the original Hamiltonian, and we will use it for the Ashkin Teller model in section 
[ill C\ Perhaps it can be useful in general, in cases where the VJs can have more than two 
possible energies. This type of operation has been used for example inlli. 
A particular case of this operation is obtained if /i = 1, and we set 
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C[ = ^Hp,) - El. (15) 



Substitution into ([131) yields 



and from (P^, the modified interaction takes the form 



^, fO ifVi(i) = i5. 

I OO otherwise 

This operation assigns infinite energy to any configuration u for which Vi{u) ^ Ei. That is, 
in the ensuing simulation the interaction V/ is frozen at energy Ei . Freezing of V/ is assigned 
with probability pi, and only when Vi{u) = Ei. It is easy to see that the SW freezing is 
precisely of this form. 



C. Cluster Method for the Ashkin Teller model 

We describe now the cluster algorithm we devised for the gAT model and the considera- 
tions that lead us to it. It consists of a freeze-delete scheme which generates non-interacting 
clusters of aligned a spins and of aligned r spins. We will phrase it in terms of the general 



scheme described in sec. IIIB 



The first decision one needs to make when coming to design a freeze-delete scheme is the 
choice of the basic interaction term V/. Our choice is to associate all interactions that reside 
on an edge < jk > of the lattice to one Vi : 

Vi = - [K^akCTj + KrTkTj + LakTkajTj]. (18) 

Since the model is invariant under any permutation of K^^, and L ( to make this symmetry 
explicit, define a new Ising spin Sj = ajTj along with the constraint sjajTj = 1 ), it is possible 
to choose L < Kcr,Kr (the reason for this choice will become clear later). The interaction 
Vi depends on four independent Ising spins that can have 16 states. Every V/ can take one 
of four possible energy values: 

El = -K„ -Kr-L 
E2 = K^-Kr + L 
E3 = -K„ + Kr + L- 
E4 = + Kr - L 

Every energy is four-fold degenerate; we denote by Ui all (four) states for which Vi = Ei. Let 
Ui represent the state of the spins aj, ak, tj, tu of a particular pair of n.n. sites < j,k >. 
Four representative possible states are depicted in figure |[ ui is a ground state in which 
both the a and the r bonds are satisfied. U2 and M3 are excited states in which one Ising 
bond is broken and the other is satisfied. Our choice L < K^, Kj. makes M4 the highest 
energy statelii, in which both Ising bonds are broken. 
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The general philosophy of our freeze-delete scheme is as follows. We wish to build clusters 
of cr spins and of r spins since we know that these clusters are the basic excitations of the 
model. Clusters of s = ar spins are not important because L, the bond between them, is 
the weakest. In order to build clusters, one needs to freeze parallel a spins to each other 
and parallel r spins to each other, and delete the bonds between antiparallel spins. For 
example in the state U2 we wish to freeze the bond between the r spins and delete the bond 
between the a spins. This consideration leads us to include in our scheme two operations. 
The operation which we will identify as i=2 freezes the bond between tj and and deletes 
the bond between aj and ak (see fig. ^ ). According to the discussion above, we want to 
perform this operation with some probability P2 ^ for the state U2, and with probability 
for the states M3 and U4. This operation allows one to move from U2 to ui, so in order to 
maintain detailed balance, we must perform it with some non vanishing probability q2 ^ 
on Ui too. All this is achieved by assigning, in the modified interaction term V2, infinite 
energy to u^jU^ and zero energy to Ui,U2- This modification is precisely of the form of eq. 
(|l^) so it is of the generalized freezing type described at the end of section pil B| and hence 
the probabilities must be (see eq. (p^)EJ: 

piw=jf"'"''i^!"i=^^°';^<"!=^'. (20) 

^ [0 It Vi{u) = or Vi{u) = E4 ^ ' 

The choice of C2 will set the values of •p2 and q2 '■ 

P2 = e^'("2)+C2 ^ ^vdui)+C2 _ (21) 

The operation which we identify as i=3 freezes the bond between aj and ak and deletes 
the bond between tj and r^. It follows the same logic as the operation i=2 and its probability 
is: 

pi(,A _ / e^'^"^^''^ if Vii^) = or Vi{u) = E, 
3^""^ ~ \ if Vi{u) = E2 or Vi{u) = ^ ^ 

Again we define the probabilities: 

^ ^vdus)+c, ^ eV,(«i)+C3 _ (23) 

The third operation we wish to perform is a freezing operation of the form defined in 
section [III B| equations (|I^-|T^). In our case it freezes the bond between aj and ak and the 



bond between tj and Tk, so it also builds the r and a clusters. We'll denote it by i=l and 
its probability is: 

^li'^) { Q otherwise ' ^"^^^ 



The fourth and last operation is the deletion operation defined in section [III B| equations 
(0,0) with the probability: 

P^{u) = e^'(")+^'*. (25) 
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These four operations define completely our freeze-delete scheme, summarized by fig. 
^. In the first column, in each row, one or two representative configurations of Ui ap- 
pear, depicting the state of aj,Tj,ak and r^. In the first row the four modified interactions 

are depicted. An upper double line for example, denotes a frozen bond between and 
(Tfc, a blank denotes a deleted bond. We have yet to determine the constants C2, C3,pi, C^. 
In the case of U4, we want to delete the bond between the r spins and the bond between the 
a spins with probability 1, as we never want an unsatisfied bond to be frozen, so we choose 
Cd = —£'4, to get Pd{ui) = 1. Having chosen Cd, all deletion probabilities are set. The rest 
of the constants are determined by the normalization condition (|^), i.e. 

= pfi ■ (26) 



Consequently, equations (pTl , P3| ) determine the constants C2, C3 and the probabilities ^2, Q's- 



Finally, pi is determined again by the normalization condition 

Pi = 1 - Pd{ui) - g2 - qs- (27) 

For completeness we'll list all the probabilities of our scheme which follow from equations 

mm-- 



Pd{u,) = e^'-^^ Vi 

p. = 1 - e^'--^* i = 2, 3 

qi = PiC^'-^' i = 2, 3 ■ 
Pi = 1 - Pdiui) -q2-qz 



(28) 



Checking that all the probabilities of the scheme fulfill the condition < Pi{u) < 1 for all 
i and u is trivial for all probabilities except for pi which is a bit more laborious, but still 
straightforward. 

We've described how we generate a new Hamiltonian of non-interacting clusters of a 
spins and of r spins. Now we have to choose some legitimate MC procedure to simulate 
this Hamiltonian. We chose to do it in a similar fashion to the single cluster algorithm of 
WolflS, but the SW version is equally applicable. To be more explicit, we choose a site j of 
the lattice at random and a random spin, either aj or Tj ( The choice of a or r was done 
with probability 1/2 which is sensible in the case K„ = K^, but in general any probability 
is acceptable, and an optimal choice can be made to minimize the autocorrelation time r). 
This spin will belong to the cluster we will flip, and the cluster will contain only r (ex) 
spins if we initially chose Tj {(Tj). We perform the freeze-delete operations only on bonds 
belonging to the surface of the cluster. For example suppose at some stage Tj was joined 
to the cluster, and suppose that the term Vj^k was modified according to i=3, then t^ will 
not be joined to the cluster and we need not perform freeze-delete operations on the other 
three bonds connecting the site k. If it is modified according to i=l or i=2 then is joined 
to the cluster. Since our cluster algorithm fits into the general scheme, we do not need to 
prove detailed balance, while ergodicity is ensured by the deletion operation. 

We've explained the reasoning behind our algorithm. We believe that the clusters we 
build of parallel a spins and of parallel r spins are the basic excitations of the model. We 
therefore believe that it will be efficient. It is encouraging to notice that in the decoupled 
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Ising subspace and the 4-state Potts subspace our freeze/delete scheme is identical with SW's 
freeze/delete operations for these models. In the appendix we show that our algorithm is 
equivalent to the idea of embedding into the AT model an Ising model and simulating it by 



Wolff's single cluster procedure. In sec. |^ we present numerical evidence for our algorithm's 
efficiency. 

D. 'Naive' SW option 

We find it illuminating to compare our algorithm to a cluster algorithm which one could 
regard as the naive generalization of the SW method to the AT model. Such an algorithm 
would define Vi in the same manner as our scheme does ( see eq. |18D. For each Ui, the 
bonds between the two neighboring sites get either deleted ( with our Pd{ui) ) or frozen with 
Pf{ui) = 1 — Pd{ui). Fig. ^ can clarify how this scheme fits into the general one of sec. 
Ill B| and how it compares with our scheme. Since it fits into the general scheme we do not 



need to prove detailed balance, while the deletion operation ensures its ergodicity except for 
T = 0. 

This scheme also generates clusters of a spins and clusters of r spins, but with the 'naive' 
scheme, antiparallel Ising spins can be found in the same cluster, so this scheme does not 
identify the elementary excitations of the model. Moreover the a clusters' structure is forced 
to match the r clusters ( and vice versa). This is totally unphysical since in practice it is 
energetically favorable for the two cluster structures not to match ( for L < K„,Kr ). At 
T = this scheme could freeze the whole lattice into a single cluster even when it is not 
in the ground state. This is in contrast to the ability of our scheme to relax excitations on 
any scale even at T = ( in other words. At T = our algorithm will bring any initial 
configuration to the ground-state in a short time). At a finite temperature the 'naive' SW 
scheme will produce clusters that are too large. In sec. |IV7\] we list results of simulations 



using the single cluster (IC) version of this 'naive' SW scheme, which show that it is indeed 
much less efficient than our algorithm. 

IV. SIMULATIONS OF THE AT CRITICAL LINE 

In addition to checking the efficiency of our method at the AT critical line, we also wanted 
to check a prediction made by Li and Sokali, who have proved a rigorous lower bound 

zsw > a/u (29) 

for the dynamical critical exponent of the SW algorithm for the ferromagnetic q-state Potts 



model. More precisely, the Li-Sokal bound is TexfJint > const x Ch, from which eq. (|29D 
follows (for definitions of Texp/int and Zexp/int seeli3). This bound relates dynamics to the 
static properties of a model and is thus of great importance. Their proof is for both Zint 
and Zexp- We wanted to check whether this bound holds all along the AT critical line, which 
connects the decoupled Ising critical point with ^ = at one end to the 4-state Potts critical 
point with ^ = 1 at the other end. At both of these points our algorithm is identical to 
Wolff's IC version of the SW method for Potts models. The Li-Sokal bound was proved for 
the SW dynamics, but perhaps it is valid for Wolff's IC dynamics as well (at least for d = 2). 
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We are unaware of any rigorous proof for that, but our results for z for q = 2 and 4 seem to 
indicate Zjn = zsw and are in accordance with previous results for the Ising critical point 
(seei and@). Besides estimating the dynamical exponents, we've estimated, using finite size 
scaling, the critical exponents ^ and ^ along the AT critical line. As stated in sec. |I| the 
AT critical line is given by 

Z =l-2X. (30) 

Our measurements were done at the decoupled Ising critical point Xq and at the 4-state Potts 
critical point X4 (see figure and exact definition in sec. ||. Three additional measurements 
were carried out at three intermediate equidistant points in the X — Z plane, on the AT 
critical line. So a total of 5 measurements were done at the points 

X, = Xo + |(X4-Xo) « = 0...4. (31) 

All the points Xj are marked in fig. [l| . We simulated lattices with periodic boundary 
conditions of up to size 128 x 128, and up to 5 x 10^ clusters were flipped for each lattice 
size. 

We calculated the energyt^ per site: 

E= {E) = --^( [KcTkCTj + KTkTj + LakTtajTj]) (32) 

^ <k,j> 

where L is the linear lattice size, and the angular brackets denote the usual thermal MC 
average. The specific heat per site follows from the energy fluctuations 

C=^ = V{{E^) - {Ef). (33) 
Kb 



We calculated the magnetic susceptibility of the cr spins defined asli 

1 

L2 



k 



and also measured c, the size of the cluster flipped at each step, and calculated (c). There 
is a connection between the size of the clusters and the susceptibilityS which is common 
for algorithms which build non-interacting clusters of spins (if all spins in a cluster have the 
same value). For IC algorithms it has the simple formtl 

X=(c). (35) 

To get the dynamic properties we calculated the time dependent autocorrelation functions 
(pEif) and 0x(t) defined as 

, _ (A(O)A(t)) - {Af 

where A stands for E or x- A typical plot of 0x(^)' measured from the X2 model at L = 128 
is presented in fig. ^ 
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A. MC results and discussion 



1. Susceptibility and Specific heat 

According to finite-size scaling theory^ one expects: 

X ~ C ~ (37) 

for large enough L. Fitting our measurements of x to eq. (|37[) for lattice sizes L > 16 fits 
the exact universal value ^ = | within errors (see table |) as expected. Our measurements 
confirm the equality (p5|), where (x) and (c) have an error of the same magnitude. Plots of 
log C vs. log L for the five models can be seen in fig. |^. From linear fits to the log — log 
plots, estimates for - are obtained which do not agree with the exact known values. For 
comparison with the exact values see table H and fig. 0. The differences may be due to 
finite size effects and corrections to scaling . An exception to this mismatch is the decoupled 
Ising point for which the value of ^ = fits nicely according to the semi-log plot in fig. 
^ (for completeness we also quote in table | an estimate for ^ at Xq from the log — log 
fit of fig. The slope of the C curves at Xi and X2 does tend to a lower value with 
increasing lattice size. The 4-state Potts model is known to have a logarithmic correction 
to scaling: C ~ L/log^''^Ll; in any case, our result - = .747(3) agrees with previous MC 
results obtained from lattices of sizes up to L = 25611. 



^exp and Tint 

In order to check the efficiency of our algorithm and to see whether the Li-Sokal bound 
holds, we measured Texp and Tint for both the energy E and the susceptibility x- This was 
done by the following procedure. (pE/xi't) was plotted on a semi- log plot, the unit of time 
being a single cluster flip. Texp was then consistently estimated by the slope of In (f){t) in a 
window from T^xp to STexp- A typical example for (f){t) and the linear fit to extract T^xp can be 
seen in figure ^ To calculate Tint, we integrated numerically (f){t) in the interval from t = 
up to t ~ l.bTexp and estimated the tail (t > Lbr^xp ) of (f){t) by the exponential fit. The 
error of both r's was estimated from repeated experiments when the statistics were large 
enough and otherwise from subjective estimates from the fluctuations in T^xp in the window 
Texp — 3re^p. Each MC step, or cluster flip, involves on the average the flipping of (c) spins. 
Since the natural unit of time is a sweep over all spins of the lattice, we multiplied r by 
(c)L-2. 

In figures |^ through |l2| we plot measured values of T^xp and Tint for E and x as a function 



of the linear lattice size L. In figures we display the measurements of at different 
models X^ separately, while in figures ^ and |T2| we display T^xp^E and Tint,E for the five 
models jointly. Fitting our results to the forms 

Texp ~ L^-'' ; Tint ~ L^*"' , (38) 

for large enough L, yielded the dynamical critical exponents listed in table |I[ Figures 
are both from the decoupled Ising model point Xq. In table | we list values of z obtained 
from the power law fits of fig. 0. None the less, comparing these fits to the logL fits of 
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fig. D, we feel that the latter ones are better and that the true value is z = for all four 
autocorrelation times r. At Xq ( and all along the Z = or L = line ) our algorithm 
is identical with the Wolff IC algorithm for the Ising model. Our results for Zint from the 
fits in fig. 1^ are consistent with those of Wolfilli. A logL behaviour has been measured for 
Texp,M and Tint,M for SW dynamics inil . 

At the point X4 , the 4-state Potts critical point (and all along the X = Z line), our 
algorithm is again identical with the Wolff IC algorithm for the 4-state Potts model. Our 
result for Zi„t b, obtained from the fit in figure O is consistent with the result of ref. i for SW 
dynamics. This result can be added to the accumulating evidence indicating: zsw = ^Woiff 
for d = 2. 



The continuous variation of z along the AT critical line can be explicitly seen in figures ^ 



and |T2| where we plot Texp,E and Tint,E for the 5 points Xj. In fig. |T2| we also plot results from 
simulations at X2 using the Metropolis method and the 'naive' SW method described in sec. 
ill L)|. For the Metropolis method, we measured a value of Zint,E = 1.64(8) , which should be 



compared with Zint,E = .542(8) using our method. The 'naive' SW method almost froze the 
whole lattice, into a single cluster, the size of which was almost independent of the lattice 
size L. For example, the average cluster size for L = 32 was ^ = .868(2) , as compared with 
|j = .481(1) using our algorithm. We found determination of r for the 'naive' SW method 
of lattice sizes L > 32 so time consuming that it was impractical. The advantage of our 
method is clearly demonstrated. The three methods yielded the same results for the static 
observables (within errors). 

The main question we wish to answer is whether the Li-Sokal bound is fulfilled and 
whether it is sharp. In fig. 13 we compare, for the 5 points Xi i = . . .4, the exact values 
of ^ and our estimated values of ^, Zint,x and Zexp,x- We see that the rise of our estimated 
values of z follows that of ^, from the decoupled Ising point Xq to the 4-state Potts point 
X4. Except for the point X4 the Li-Sokal bound is fulfilled with respect to {-)exact- The 
anomalously low values for z at X4 are probably caused by the multiplicative logarithmic 
correction for C described in the discussion of the specific heat results. This explanation 
has been suggested by Li and Sokali to explain the low value of zsw for this model. From 
table H we see that z > estimate with only one exception at the point Xi (see fig. ^. 
The fact that Zint,x < estimate could be due to a difference in the finite size corrections 
or in the corrections to scaling. Tint,x is probably less infiuenced by one or both of these 
two factors, in comparison with the large difference between {^)estimate and {^)exact- At the 
point X3, where the difference between {^)estimate and (^)exact is the smallest, no anomahes 
in the Li-Sokal bound occur (see fig. [Tol) . 

We conclude that the Li-Sokal bound is fulfilled in a moderately sharp manner. The 
smallest z (which is Zint,x) fulfills {^)estimate < Zint,x < {j;)estimate + 0.1 (uot including the 
anomalies of the models X4 and Xi discussed above) . Note that we are comparing estimated 
values of z with estimated values of - . This is the correct comparison to make, assuming 
that r and C have similar finite size corrections and similar corrections to scaling. 



V. SUMMARY AND DISCUSSION 

The correct identification of the basic excitations of the model lead, along with the guide- 
lines of a general scheme for cluster algorithms, to the construction of a cluster algorithm for 
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the AT model. The algorithm was shown to be identical with the one obtained by embed- 
ding Ising spins into the AT model. Our algorithm is suitable for spatially varying coupling 
constants under the restriction L'''^ < K^'\ K^'^ everywhere on the lattice. We are currently 
carrying out intensive simulations of a random-bond version of the AT model. 

The dynamical behaviour of the cluster algorithm was examined on the AT critical line. 
Critical slowing down of the algorithm ( < 2; < 1. ) was found to be significantly smaller 
than that of the standard Metropolis method. A continuous variation of the dynamical 
exponent z along the AT critical line was seen, along with continuous variation of the static 
exponent -. The Li-Sokal bound 2 > -, that was proved only for g-state Potts models with 
SW dynamics, is satisfied for the AT model with single cluster dynamics. The bound is 
moderately sharp. Another new result is zwoiff = zsw for the 4-state Potts model, which 
can be added to the accumulating results indicating that in two dimensions zwoiff = zsw ■ 
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APPENDIX: ISING EMBEDDING 

Our algorithm can also be seen from a totally different point of view; as an example of 
an embedding algorithmi'0. The main idea is to embed into the AT model an Ising model 
of space dependent couplings Jjk and simulate it using the SW or Wolff procedure for the 
Ising model. To be explicit, consider the Hamiltonian (|), and take the r variables as fixed, 
so we can write 

H = Hl + H2 = - i^o- + LTkTj)akCrj - KrTkTj. (Al) 

<k,j> <k,j> 



T1 



represented by the second sum is a constant, and remembering that we chose L < K„, K. 
Til is a ferromagnetic Ising model in the a variables with couplings Jjk = K„ + Lr^Tj. 
Simulating this Ising model with any procedure that will maintain detailed balance with 
respect to Tii and will not change the value of 7^2, will also maintain detailed balance 
with respect to Ti. So we can use for example the SW or Wolff procedure for the Ising 



model, explained in sec. |III B| . This by itself will maintain detailed balance but will not 
be ergodic since the r variables will not be updated. Obviously to update the r variables 
the same process should be repeated , holding the a variables fixed and simulating an Ising 
Hamiltonian with the r variables. To summarize, a possible procedure would go as follows: 
Choose at random whether to embed into the AT Hamiltonian an Ising Hamiltonian in the 
a spins or in the r spins. Pick a random site in the lattice, grow a cluster of cr or r spins 
using the Wolff (IC) procedure using the Ising Hamiltonian and flip it. As we will now 
show, this process ( we'll call it for shortness IE - Ising embedding) is exactly equivalent to 
the AT algorithm (ATA) we suggested in section |III C . 



A first argument which is really sufficient goes as follows. Suppose we perform the IE 
freeze-delete scheme for the a spins for the whole lattice, viewing the r spins as constant. We 
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could say we've generated a new Hamiltonian of noninteracting clusters of a spins and which 
assigns infinite energy to any configuration which differs from the current configuration of 
T spins. Now suppose we perform the ATA freeze-delete scheme for the whole lattice. Then 
we get a new Hamiltonian of noninteracting clusters of a spins and of r spins, but if at 
this stage we decide to flip only a spins, then in practice we've assigned infinite energy to 
any configuration which differs from the current configuration of r spins. In practice we flip 
only one cluster of the a spins, so what we do is identical to the Wolff (IC) version of the 
IE. Since both processes maintain detailed balance using the same new Hamiltonian they 
must do so using the same probabilities (the general scheme in sec. |111B shows a one to 



one correspondence between probabilities and the new Hamiltonian.), which completes our 
argument that the two methods are actually identical. 

One can, of course, check that the probabilities of the two procedures are the same. For 
example, denote the probability to delete a bond between a spins in ATA as P/J^(-u). for 
Ui for example. From fig. ^: 

i^^r («i) = Uui) + q2 = e^-^^ = e-2(^^+^). (A2) 
Now with the Ising Embedding algorithm, according to (^: 

Pl^iu) = e-(^-+i-*-.)(-fe-.+i), (A3) 



so Pd^{ui) = P^J (ui). This check can be carried out for all 



Ui. 
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FIGURES 



FIG. 1. Phase diagram of the Ashkin Teller (Z(4) ) model. 

FIG. 2. States of spins at a pair of nearest-neighbor sites j, A;. Each state Ui represents one out 
of four states with the same energy £'j. 

FIG. 3. Freeze/delete scheme for the AT model. The four states of two neighboring sites 
are denoted in the left column; Four freeze/delete operations are denoted in the top row, and the 
respective probabilities arc denoted in the table. A frozen bond is denoted by a double thick line, 
while a deleted bond is denoted by a blank. 

FIG. 4. 'Naive' SW freeze/delete scheme for the AT model. The four states of two neighboring 
sites are denoted in the left column; Four freeze/delete operations are denoted in the top row, and 
the respective probabilities are denoted in the table. A frozen bond is denoted by two thick lines, 
while a deleted bond is denoted by a blank. 

FIG. 5. A log-log plot of specific heat vs. linear lattice size at the five critical models Xq . . . X4. 
The solid lines are the linear fits, the critical exponents ^ are given in table 1. 

FIG. 6. The In of the time auto-correlation function of the susceptibility (py-{t), measured for 
the X2 model at L=128. The unit of time is a single MC step or a single cluster flip. The vertical 
lines are the error bars. The exponential fit is also shown. 

FIG. 7. Log-log plots of auto-correlation times and the specific heat C as a function of 
linear lattice size L, at the decoupled Ising point Xq. The values of z listed in table 1 are the slopes 
of the fits. 

FIG. 8. Semi-log plots of auto-corrclation times and the specific heat C as a function of 
linear lattice size L, at the decoupled Ising point Xq. This fit seems to be better than the log-log 
fit, yielding values of z = ^ = 0. 

FIG. 9. Log-log plots of auto-correlation times and the specific heat C as a function of 
linear lattice size L, at the point Xi. The values of z listed in table 1 are the slopes of the fits. 

FIG. 10. Log-log plots of auto-correlation times Ty. and the specific heat C as a function of 
linear lattice size L, at the point X3. The values of z listed in table 1 are the slopes of the fits. 

FIG. 11. A log-log plot of Texp,E VS. linear lattice size at the five critical models Xq . . .X4. 
The solid lines are the linear fits, the critical exponents Zexp,E are given in table 1. 
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FIG. 12. Log-log plots of Tint,E as a function of linear lattice size L, at all 5 points i = . . . 4. 
The continuous variation of the slope z along the AT critical line is easily seen. Metropolis results 
and 'naive' SW results for X2 are also plotted. 

FIG. 13. Gomparison of z and ^ at the 5 points Xi i = . . .4. The line denotes the exact 
value of ^, while our estimated values for the five models are denoted by full circles, values of 
Zint,x/ ■^exp,x denoted by empty circles/ crosses. 
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TABLES 



TABLE I. Results from the AT critical line. The errors in parentheses are only the statistical 
errors of the fit in our fitting interval. They do not include systematic errors stemming from finite 
size effects and corrections to scaling. 





1 


a 

V 










a 


^0 


1.751(4) 


.23(1) /O 


.26(3) 


.23(2) 


.13(1) 


.267(4) 





Xi 


1.751(1) 


.38(2) 


.396(5) 


.40(3) 


.273(3) 


.47(2) 


.2096 


X2 


1.752(3) 


.542(8) 


.61(1) 


.68(3) 


.532(5) 


.62(1) 


.4182 


X3 


1.763(6) 


.655(13) 


.719(2) 


.74(2) 


.717(8) 


.72(3) 


.6383 


Xi 


1.752(4) 


.747(3) 


.92(1) 


.99(5) 


.99(3) 


.94(3) 


1 
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